!> author: 左志华
!> date: 2022-06-17
!> license: BSD 3-Clause
!>
!> Density summation <br>
!> 针对实粒子进行密度求和：<br>
!> 1. 连续密度法；<br>
!> 2. 正则密度求和法；<br>
!> 3. 密度求和法。
module sph_density_summation

    use sph_kinds, only: rk
    use sph_particle, only: particle_t
    use twin_sph_list_type, only: twin_list_t
    use sph_smoothed_kernel_function, only: skf

    integer :: density_summation_type   !! density summation type <br>
                                        !! 密度求和类型：<br>
                                        !! 1. 连续密度法；<br>
                                        !! 2. 正则密度求和法；<br>
                                        !! 3. 非正则密度求和法
    procedure(density_summation_fcn), pointer :: density => continuous_density_summation
                                                                !! Density summation function pointer <br>
                                                                !! 密度求和函数指针

    abstract interface
        !> Density summation <br>
        !> 密度求和
        subroutine density_summation_fcn(particle, pairs, n)
            import :: particle_t, twin_list_t
            type(particle_t), intent(inout) :: particle             !! particle <br>
                                                                !! 区域
            type(twin_list_t), intent(in) :: pairs              !! particle pairs <br>
                                                                !! 粒子对
            integer, intent(in) :: n                            !! ntotal <br>
                                                                !! 总粒子数
                                                                !! @note 用以初始化本地数组长度
        end subroutine density_summation_fcn

    end interface

contains

    !> Density summation <br>
    !> 对所有粒子进行非正则化密度求和 <br>
    !> @note 更改内容：实虚粒子的密度
    pure subroutine denormalized_density_summation(particle, pairs, n)
        type(particle_t), intent(inout) :: particle             !! particle <br>
                                                            !! 区域
        type(twin_list_t), intent(in) :: pairs              !! particle pairs <br>
                                                            !! 粒子对
        integer, intent(in) :: n                            !! ntotal <br>
                                                            !! 总粒子数 @note 未实际使用
        integer :: i

        particle%rho(:) = skf%selfden*particle%mass(:)              ! 计算自身对密度的影响，并初始化当前步长密度值

        do concurrent(i=1:pairs%log)                        ! 计算密度求和
            associate (ip => pairs%items(1, i), jp => pairs%items(2, i))

                particle%rho(ip) = particle%rho(ip) + particle%mass(jp)*pairs%w(i)
                particle%rho(jp) = particle%rho(jp) + particle%mass(ip)*pairs%w(i)

            end associate
        end do

    end subroutine denormalized_density_summation

    !> Normalized density summation <br>
    !> 对所有粒子进行正则化密度求和法
    !> @note 计算量较密度求和法大
    pure subroutine normalized_density_summation(particle, pairs, n)
        type(particle_t), intent(inout) :: particle !! particle <br>
                                                !! 区域
        type(twin_list_t), intent(in) :: pairs  !! particle pairs <br>
                                                !! 粒子对
        integer, intent(in) :: n                !! ntotal <br>
                                                !! 总粒子数
                                                !! @note 用以定义本地数组长度
        integer :: i
        real(rk) :: rho_tmp(n), wi(n)           !! 分母部分

        rho_tmp(:) = particle%rho(:)              ! 临时存储密度值

        particle%rho(:) = skf%selfden*particle%mass(:)
        wi(:) = particle%rho(:)/rho_tmp(:)

        do concurrent(i=1:pairs%log)            ! 计算密度求和
            associate (ip => pairs%items(1, i), jp => pairs%items(2, i))

                particle%rho(ip) = particle%rho(ip) + particle%mass(jp)*pairs%w(i)
                particle%rho(jp) = particle%rho(jp) + particle%mass(ip)*pairs%w(i)
                wi(ip) = wi(ip) + particle%mass(jp)*pairs%w(i)/rho_tmp(jp)
                wi(jp) = wi(jp) + particle%mass(ip)*pairs%w(i)/rho_tmp(ip)

            end associate
        end do

        particle%rho(:) = particle%rho(:)/wi(:)  ! 正则化密度

    end subroutine normalized_density_summation

    !> Continuous density summation <br>
    !> 基于连续性方程的连续密度法，计算密度的增量以在时间积分算法中更新密度值
    !> @todo itype 相同才进行密度求和？
    pure subroutine continuous_density_summation(particle, pairs, n)
        type(particle_t), intent(inout) :: particle !! particle <br>
                                                !! 区域
        type(twin_list_t), intent(in) :: pairs  !! particle pairs <br>
                                                !! 粒子对
        integer, intent(in) :: n                !! ntotal <br>
                                                !! 总粒子数
                                                !! @note 未实际使用
        integer :: i
        real(rk) :: dvel(2)

        particle%drho(:) = 0.0_rk
        do concurrent(i=1:pairs%log)
            associate (ip => pairs%items(1, i), jp => pairs%items(2, i))

                dvel(:) = particle%vel(:, ip) - particle%vel(:, jp)

                associate (vcc => sum(dvel(:)*pairs%dwdx(:, i)))
                    particle%drho(ip) = particle%drho(ip) + particle%mass(jp)*vcc
                    particle%drho(jp) = particle%drho(jp) + particle%mass(ip)*vcc
                end associate

            end associate
        end do

    end subroutine continuous_density_summation

end module sph_density_summation
